using System;
using System.Reflection;
using NBody;
using NBody.Analysis;
using OptionParser;
using DataFile;

[assembly: AssemblyTitle ("NBodyPositions")]
[assembly: AssemblyVersion ("1.0.*")]
[assembly: AssemblyDescription ("Output the positions of particles in an N-body system.")]
[assembly: AssemblyCopyright ("2005 Joseph D. MacMillan")]

public class PositionsOptions : InputOutputGetOptions
{
    [Option ("-c", "Adjust for c/m")]
    public bool CentreOfMass;
    
    [Option ("-x", "Rotate about x axis (degrees)")]
    public double XAngle;
    
    [Option ("-y", "Rotate about y axis (degrees)")]
    public double YAngle;
    
    [Option ("-z", "Rotate about z axis (degrees)")]
    public double ZAngle;
    
    [Option ("-m", "Rotate by matrix in file")]
    public string MatrixFile;
    
    [Option ("-n", "Output n positions")]
    public int N;
    
    [Option ("-r", "Output particles less than radius")]
    public double R;
    
    public PositionsOptions()
    {
        XAngle = 0.0;
        YAngle = 0.0;
        ZAngle = 0.0;
        MatrixFile = "";
        N = 50000;
        R = Double.PositiveInfinity;
    }
}

public class NBodyPositions
{
    public static void Main(string[] args)
    {
        PositionsOptions opts = new PositionsOptions();
        opts.ProcessArgs(args, ".pos");
                        
        NBodySystem s = NBodySystem.Read(opts.InFile);
        
        if (opts.CentreOfMass)
            s.AdjustForCentreOfMass(1e-3);
        
        if (opts.R < Double.PositiveInfinity)
            s.ExtractSphere(opts.R);
        
        if (opts.XAngle != 0.0)
            s.RotateAboutX(opts.XAngle);
        if (opts.YAngle != 0.0)
            s.RotateAboutY(opts.YAngle);
        if (opts.ZAngle != 0.0)
            s.RotateAboutZ(opts.ZAngle);
            
        if (opts.MatrixFile != "")
        {
            Table a =  Table.Read(opts.MatrixFile);
            s.RotateByMatrix(a.Array);
        }
        
        double mod = (int)Math.Ceiling((double)s.NumParts / (double)opts.N);
        opts.N = (int)Math.Ceiling((double)s.NumParts / (double)mod);
        int mod_value = (int)Math.Ceiling((double)s.NumParts / (double)opts.N);
        
        Table data = new Table(opts.N, 3, "Particle positions (x, y, z)");
        int c = 0;
        for (int i = 0; i < s.NumParts; i++)
        {
            if (i % mod_value == 0)
            {
                data[c, 0] = s[i].X;
                data[c, 1] = s[i].Y;
                data[c, 2] = s[i].Z;
                c++;
            }
        }               
        
        data.Print(opts.OutFile);
    }
}
